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Jj~H^ We present a molecular dynamics and kinetic theory study of granular material, modeled 

_"Ji ■ by inelastic hard disks, fluidized by a random driving force. The focus is on collisional averages 

and short distance correlations in the non-equilibrium steady state, in order to analyze in a 

C*^ , quantitative manner the breakdown of molecular chaos, i.e. factorization of the two-particle 

' distribution function, /' '{x\,X2) — xf^){^i)f (2^2) in a product of single particle ones, 

(-H ' where Xi — {ri,Vi} with i = 1, 2 and x represents the position correlation. We have found 

fi . that molecular chaos is only violated in a small region of the two-particle phase space {ii, X2}, 

Q^ ' where there is a predominance of grazing collisions. The size of this singular region grows with 

C [ increasing inelasticity. The existence of particle- and noise-induced recoUisions magnifies the 

I ■ departure from mean field behavior. The implications of this breakdown in several physical 

"t^ ' quantities are explored. 
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I. INTRODUCTION 



The interesting phenomena observed in recent experiments with mono- and multi-layers of granular material 
' ^ on vibrating plates ||l|-y show the need to develop kinetic theories for rapid granular flows with mechanisms 
for energy input, different from those in shear flows or flows through vertical pipes. In the present article the 
fluidization is driven by a random external force, which gives frequent kicks to each particle in between collisions. 
Such a driving mechanism has recently been studied by many authors p|-pT[| . The basic physical interest is the 
. understanding the non-equilibrium stationary states (NESS) that exist in the presence of this random force. 
^s^J ' The advantage of this fluidization mechanism, besides its potential physical realizations, lies in the fact that 
J^ the NESS is linearly stable against spatial inhomogeneities. 



r — In Ref. 1 10 1, to which we refer as paper I, we have studied the large scale structure and presented a hydro- 

l/~j ' dynamic description of randomly driven granular fluids, modeled as systems of smooth inelastic hard spheres 

C — ' (IHS). The IHS model accounts for two essential features of granular matter: hard core exclusion and dissipa- 

^D ' tive collisions II4] . The dynamics is described by a constant coefficient a of normal restitution. In collisions a 

J^ fraction of the relative kinetic energy is lost, which is proportional to the inelasticity e — \ — o? . The stochastic 

^^ ' external force compensates this energy loss, and drives the IHS fluid into a NESS. This stationary state, though 

homogeneous and stable against spatial fluctuations on large space and time scales (at least for weakly inelastic 

spheres), was shown to exhibit long range spatial correlations in density, velocity and granular temperature 

, fields that extend much beyond the mean free path. In fact, the corresponding structure functions S'(A;) diverge 

i-rt ' as 1/A:^ as the wave number fc — > 0, a behavior caused by the random external force, which does not conserve 

P5 momentum whereas the collisions between particles do. These long range correlations are of algebraic form, 

O ^ l/r"^^^, which corresponds to Inr in two dimensions (d — 2). The existence of such extremely long range 

J-J spatial correlations is one example of the many nontrivial properties of non-equilibrium stationary states in 

> ; general ^Q. 

l^ • Differences in the stationary states between fluids with dissipative and conservative interactions also manifest 

rS themselves in the kinetic properties of the fluid, such as the velocity distribution function, which deviates 

j^ from a Maxwellian in particular in the high energy tail of the distribution. In Ref. p| the existence of an 
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overpopulated high energy tail / ~ exp[—Cv'^^'^], where C is a constant that depends on the inelasticity, has 
been obtained from kinetic theory. A similar behavior has been observed experimentally at high vibrational 
accelerations ^g- This observation indicates that certain features of the experiment might be reproduced 
by modeling the input of energy into the horizontal motion of the beads by a random external force, although 
other energy injection mechanisms that could be relevant to recover the large velocity tail have been put forward 
|15| . In similar experiments with a vertically vibrating plate covered with a mono-layer of steel balls with a 
packing fraction around 50% the velocity distribution of the horizontal velocities has been measured, and again 
overpopulated non-Gaussian high energy tails have been observed. In the present paper we will investigate the 
kinetic properties and short scale correlations that characterize the stationary state. More specifically, we will 
compare molecular dynamics (MD) simulations of inelastically colliding disks with analytic predictions based 
on the assumption of molecular chaos. 

The Boltzmann equation for dilute gases of particles that interact via short-ranged repulsive interactions is 
based on the assumption of molecular chaos, also called the Stosszahlansatz or mean field approximation. It 
assumes that the velocities of colliding particles just before collisions are uncorrelated, i.e. their pair distribution 
function factorizes, f^'^'{xi,X2,t) = f{xi,t)f{x2,t), where Xi = {rj,Vj} denotes the position and velocity of 
particle i. Enskog's extension of the Boltzmann equation to a dense system of hard spheres (l6|, referred 
to as Enskog-Boltzmann equation, is also based on the fundamental assumption of the absence of velocity 
correlations. Here the assumption of molecular chaos postulates that f^'^''{xi,X2,t) = xf{^iTt)f{x2,t) for 
approaching particles (vi2 ■ ri2 < 0) just before collision (ri2 = a -|- 0), where x is assumed to be the radial 
pair distribution function at contact g{ri2 = cr + 0) in local equilibrium. It implies the additional assumption 
that spatial correlations between colliding particles just before collision are independent of their velocities, 
i.e. absence of position-velocity correlations. The Enskog x factor enhances the collision frequency at higher 
densities. For dilute gases the assumption of molecular chaos seems to be justified. Recently Lutsko Q and 
Soto and Mareschal llq] derived for a freely evolving inelastic hard disk fluid a relation between pre- and post- 
collision radial distribution function at contact, as a function of the angle, 9 = cos~^(vi2 • fi2), between the 
relative velocity Vi2 of the colliding particles, and their relative position at contact ri2, and they confirmed 
their results by MD simulations. Their observations made it clear that further arguments are needed to clarify 
the meaning of the %— factor in Enskog's formulation of the molecular chaos assumption. This will be done in 
section II A. 

The breakdown of molecular chaos at higher densities in classical fluids with conservative forces has been 
extensively investigated in the sixties and seventies |jl9] . This breakdown is caused by sequences of correlated 
binary collisions, the so called ring collisions EQ]. They lead to long time tails in velocity and stress autocorre- 



lation functions | |2l| , p2| , and to long range spatial correlations in NESS |13|. The quantitative effects of velocity 
correlations on transport coefficients at liquid densities are also significant. For instance, molecular dynamics 
simulations on elastic hard sphere systems at liquid densities [p3] have shown that the long time tails increase 
the measured self-diffusion coefficient D typically by 15% to 20% with respect to the prediction of the Enskog 
theory, De = -Db/x, where Db is the Boltzmann value of the self-diffusion coefficient. 

A well-known example of short scale structure in granular fluids are the position- velocity correlations leading 
to the phenomenon of inelastic collapse [^,^ , which is a divergence of the collision frequency oj in a finite 
time. The collapse singularity implies that an infinite number of collisions occurs within a finite time interval 
in a subset of (nearly) touching particles, ordered in linear arrays. The phenomenon is, however, an artifact of 
the assumption that the coefficient of restitution a is independent of the impact velocities, whereas on physical 
grounds a{vi2) -^ 1 (elastic limit), as the relative velocity wi2 vanishes. Molecular dynamics simulations have 
shown that the assumption of molecular chaos is also violated in undriven granular fluids in their late stages 
of evolution, the so called nonlinear clustering regime. For instance, the measured distribution of impact 
parameters is not uniform, as expected on the basis of molecular chaos, but biased toward grazing collisions 
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As shown below, in the driven IHS fluid there is an important additional reason for the breakdown of molecular 

chaos, namely the strong increase in collision frequency at small relative velocities between two isolated particles, 
caused by the so-called noise induced re-collisions. This correction to the collision frequency, that is important 
at all densities, is also neglected in the molecular chaos assumption. 

The main goal of this paper is to quantify, analyze and interpret the effects of the breakdown of molecular 
chaos in the NESS of inelastic hard spheres that are subject to a random external force between collisions. 
We will focus in particular on velocity-velocity correlations and position-velocity correlations between particles 



almost in contact, i.e. the short scale structure. 

Section || presents the analytic results, based on the Enskog-Boltzmann equation, which has been modified to 



account for the external energy input. In section [11 we present molecular dynamics results for several quantities 
that characterize the collision processes and related short scale structure of the NESS, and make a comparison 
with predictions based on molecular chaos. 

II. KINETIC THEORY FOR THE NESS 

A. Molecular chaos and Enskog approximation 

The Enskog-Boltzmann equation for the single-particle distribution /(vi,i) in a spatially homogeneous ran- 
domly driven fluid of inelastic hard spheres of diameter tr reads in d = 2 or 3 dimensions |S[| : 

5t/(vi,t) =nxfT'*"^/dv2/do-e(vi2 • (t)(vi2 ■ cr) x 

{7^/(vr,i)/(vr,t)-/(vi,i)/(v2,i)} + f (9l7)'/(vi,t), (1) 

where Vi2 = vi — V2 and n the number density. The Heaviside function Q{x) restricts the <t integration 
to the hemisphere Vi2 • o- > 0, where <t is the unit vector along the line of centers of the colliding spheres 
at contact, pointing from particle 2 to 1. In the sequel a — a/|a| denotes a unit vector. The gain term of 
the collision integral describes the restituting collisions that convert the precoUision velocities (v|*,V2*) into 
(vi, V2), while the loss term describes the direct collisions, and contains the precoUision velocities (vi, y^) leading 
to postcoUision velocities (v^ , Vj ) . The postcollision and restituting velocities have been defined in pfl] . The x 
factor will be discussed below. 

As derived in |p|, the Fokker-Planck term accounts for the external energy input, and describes diffusion 
in velocity space with a diffusivity proportional to the rate of energy input ^q per unit mass. Here ^0 is the 
strength of the random driving force, which is assumed to be Gaussian white noise P,[lO|. 

Before studying the short scale structure, we consider the single particle distribution function /(v) in the 
NESS. The stationary solution of (|l|) is determined by the balance, to^q = F, of external heating, to^q, and 
internal loss of energy due to collisions, T. It is characterized by a time independent temperature, T = {mv^ /d), 
defined as the average kinetic energy per particle, and discussed in paper I. As mentioned in the introduction, 
this stationary solution exhibits an overpopulated high energy tail / ~ exp[— Cw"^/^]. The structure of the tail 
distribution is determined by collisions of very energetic particles with 'thermal' particles, and can be obtained 
by neglecting the gain term in the Boltzmann equation ||9(]. 

In Ref. /(v) has been calculated by solving the Enskog-Boltzmann equation (|^) by an expansion in Sonine 
polynomials. To formulate this result it is convenient to introduce a rescaled distribution function /(c), defined 
by /(v) = [l/uo]/(c) with c = v/wq, where vq = \J2T jm is the thermal velocity and d the dimensionality. This 
gives 



/(c) = ^(c) 1 + a2 



ic4-i(d + 2)c2 + id(d + 2) 



(2) 



where the Maxwellian Lp(c) = tt ''/^ exp(— c^). Note that 02 is proportional to the fourth cumulant of the scaling 
form /(c), i.e. 

«2 = ^^^ [(c^) - \d{d + 2)] = I [(4) - 3(4)^] , (3) 

and vanishes in the clastic limit. An explicit calculation to linear order in 02 gives [Q 

„, ^ 16(l-^)(l~2a^) 

73 + 56rf-24arf- 105a + 30(1 -a)a2' ^' 

In the next section this prediction will be tested against molecular dynamics simulations. 



Consider first the exact expression for the mean colhsion frequency in the homogeneous NESS, defined as 

w = na'^-^ / dvi / dv2 / dCTe(-vi2 • 6-)|vi2 • 6-\P'>{vi,V2,(t), (5) 

where /'^^'(vi, V2,(t) is the dynamic or constrained pair distribution function with velocities aiming to coUide, 
just before contact with ri2 = cr. Molecular chaos, also referred to as mean field theory, requires the complete 
factorization of the dynamic precoUisional pair function, 

/(2)(vi,V2,(t) = x/(vi)/(v2). (6) 

What is the meaning of the x— factor, used in formulating the molecular chaos hypothesis? This hypothesis for 
dilute gases concerns the absence of correlations in precollision velocities, and in precollision positions (x = 1). 
In dense fluids on the other hand, the precollision position correlation x is different from 1, but the precollision 
velocity-velocity and position-velocity correlations are still assumed to be absent. 

In the literature it is common to take x equal to the radial distribution function at contact in local equilibrium, 
i.e. X = xe = 5cq('' —^ cr -f 0), which mainly accounts for precollision hard core exclusion effects. For hard disks 
and hard spheres the latter function is approximately given by the Verlet-Levesque (2D) and Carnahan-Starling 
(3D) approximations |3l| ], 

Xe(0) = (1 - ^</')/(l ~ </-)' (2D) 

Xe(0) = (2-(/-)/2(1-0)3 (313)^ (7) 

where the packing fraction in d dimensions is defined as cj) ~ n{a/2)'^Qd/d, and Q^i — 2TT'^^^/r{d/2) is the surface 
area of a d-dimensional unit sphere. In this paper we refer to the molecular chaos approximation with x = Xe, 
as the Enskog approximation. 

In principle, different options are open for the x~factor. As f^'^' is the dynamic precollision pair distribution 
function at contact, an alternative choice for the x in the factorized form could be the dynamic precollision 
radial distribution function at contact, defined as an average over the precollision hemisphere. 



[2/nd] J dvi / dv2 f d^e(-vi2 • o-)/(2) (vi, V2, a- 



X^-^ = [2/^d]j dwi j dw2 I d^e(-Vi2-o-)/^^J(vi,V2,cr). (8) 

Another option could be the unconstrained radial distribution, g{r), in the NESS, extrapolated to contact 
(r ^ cr). This function is further discussed in subsection II D. 

For the randomly heated fluid under study here, the dynamics is not purely hard-sphere like. The random 
force acting on the particles may be expected to contribute to the value of the pair distribution function at 
contact. This effect will be addressed in the subsequent sections. 

Equation (||) with /(v) replaced by the Maxwellian, yields for the collision frequency in the molecular chaos 
approximation {jj-aic{T) — xwo(r) , and more specifically in the Enskog approximation, 

^e(T) = XEC^o(T). (9) 

Here the Boltzmann collision frequency for dilute gases is given by. 



cjo(r) = nana'^-^^Tlnm, (10) 

and the small correction of C(a2) appearing in (g) has been neglected. 

Spatial correlation functions in non-equilibrium stationary states are quite different from local equilibrium 
ones, and show l ong range correlations due to correlated sequences of ring collisions, also referred to as mode 
coupling effects |l3[ . In paper I we have shown the existence of very long range correlations ^ \jr^~'^ in the 
randomly driven IHS fluid. The short range correlations in the NESS can in principle be obtained from the 
ring kinetic equation for IHS (see Ref. |2£|). However, systematic methods to evaluate collision integrals and 
pair correlation functions at short distances using this ring kinetic theory have not yet been developed. In 
the section on sinnilation results, we return to the effects of ring collisions, and present arguments why their 
contributions are expected to be more important here than for elastic hard spheres. 



B. Collisional averages 

In hard sphere systems there are many properties that mvolvc the pah distribution function of particles just 
before coUision. To study these, it is convenient to introduce the collisional average (...)coii for a quantity A in 
the NESS, defined as 

/dvi/dv2/do-e(-Vi2 •0-)|vi2 •(T|A(vi,V2,Cr)/(2)(vi,V2,(T) 

(A(Vi, V2, cr))coIl = ^^^ • (11) 

Jdvi /dV2/do-0(-Vi2 • 0-)|vi2 • 0-|/(2)(vi, V2, cr) 
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In the sequel it is more convenient to work with a rescaled pair distribution function /^^^(vi, V2, ri2) — [l/v, 
/(^^(ci, C2, ri2). To express the collisional averages (|ll|) in rescaled variables one replaces Z*^^^ by f^^\ Vi by Ci, 
V12 by g = Ci -C2, and A(vi,V2,cr) by ^(woCi, woC2,cr). 

These objects can be conveniently computed in event driven molecular dynamics algorithms for hard sphere 
systems Q. Collisional averages are defined for particles that are about to coUide (i.e. |ri2| = tr + 0), and 
can be calculated from kinetic theory using the molecular chaos assumption, possibly supplemented with the 
Enskog approximation at higher densities. 

Collisional averages of great importance are the collisional energy loss per unit time, o'^r, and the excess 
hydrostatic pressure, p — nT, resulting from collisional transfer of momentum. With a minor generalization to 
d dimensions we obtain from Ref. P4] the exact expression for the pressure in the NESS: 



p(T)_^^ fl + a 



nT \ 2d 

1 + a\ (7LU 



na" J dci J dc2 J d&Q{-g- &)\g- &\^f^^Hci,C2, a) 

•^DcoU. (12) 



2d J vo 

The second equality is obtained by introducing the collisional average ( [ll[ ) and expressing its denominator in 
terms of the collision frequency given by (H). In fact, inserting (ra) into the first line of dl2) allows one to 
carry out the a— integration, and the right hand side becomes proportional to the rescaled velocity average 
Jdci J'dc2(?^/(ci)/(c2) = 2 without any further assumption about neglecting the term proportional to 02 in 
(||). This argument is special for the pressure, as other moments involve the knowledge of the complete /(c). 
Indeed, the generic collisional average becomes (|g • (t|'")co1i = 2™/^ T{^m + 1), independent of dimensionality, 
assuming molecular chaos (g) and replacing /(c) by the Maxwellian (p{c) (the contributions coming from 02 are 
quite small and can be neglected; they have been computed in M). Finally, the pressure can be expressed as, 

^""^^^-l = 2'^-2(l + a)x0. (13) 



nT 

Different choices for x yield different approximations. For instance, with x = Xe we obtain the Enskog approx- 
imation pe{T) for the pressure of IHS. 

In the elastic limit PEiT) at a = 1 gives the standard equation of state for elastic hard disks or spheres. Notice 
that the pressure for IHS is only defined kinetically as the momentum flux, which leads to (p^. A statistical 
mechanical derivation of the equation of state from the partition function or free energy for the IHS fluid does 
not exist. 

In a similar manner we obtain the exact expression for the collisional damping rate, 

r(r)= n-^\na''-\oTjdc,Jdc2jd&Q{-g-&)\g-a\^f^'Hc,,C2,a) 

= 7owr(|g • (Tp)coii 

= meo, (14) 

where 70 = (1 — a^)/2d is the dimensionless damping constant introduced in Refs. ||9|,[l0|]. The last equality ( [l4[ ) 
expresses the balance between the energy input due to the white noise, and the collisional loss of energy in the 
NESS, and determines the temperature T in the NESS. By specializing this equation to the Enskog mean fleld 
approximation, f^'^' = xsf f, where (|g • o-p)coU = 2, we obtain the approximate result. 



Te{T) = 27oc^e(T)T, (15) 

and similar relations for different choices of x- It is convenient to define a reference temperature Te through 
the relation, 

rE(rE) = m^o, (16) 

or more explicitly 

lE — m [ - — — T— r . (17) 

Moreover the definition of Te combined with the NESS condition V{T) = to^q implies the relation V{T) = 
rE(rE), and consequently, 

r(r) _rE(rE) _ fTE\"\ ^^^^ 



Te{T) Te{T) V T 
In the sequel we will also use a reference frequency uje without any argument to denote. 



UE = uJe{Te) = XE^dna'^-WTE/TTm. (19) 

Although we will discuss the simulation results in detail in the next section, it is of interest already at this 
point to note that for these systems, as shown in Fig. n^, the ratio of the kinetic temperature and the reference 
temperature, T/Te, is only somewhat larger than 1 for all a, that it approaches 1 in the elastic limit (a — > 1), 
and that it monotonically increases with decreasing a (see Fig. M. The same figure shows that the ratio, uo/loe 
also approaches 1 for a -^ 1, with a steep increase to a value 5.6 as a — > 0. Further discussion of these points is 
postponed till section III. If the Enskog factorization, /(^) = Xe//, would be exact, then T — Te and oj — loe- 
A third quantity of interest, the precoUisional x^"'' factor, defined in (jg), can also be expressed as a coUisional 
average using (pl|), 

X^-) = (2c^/firfna'^-i«o)(|g • ^l^)coiv (20) 

Before closing this section a caveat about internal consistency is appropriate. To obtain consistent theoretical 
predictions for the pressure p or dissipation rate F, it is paramount that both factors, w and (|g • cr|'")coii, be 
calculated using identical approximations for f^"^' . For instance, the mean field or molecular chaos approximation 
for the dissipation rate, rE(r) — 27owT, - an expression commonly used in granular hydrodynamic equations 
- should necessarily be combined with loe{T) in (0). Any improved theoretical calculation for lo without a 
concomitant correction to the mean field result for (|g • cr|™)con is inconsistent . 

C. Velocity distributions 

We study a variety of coUisional averages and corresponding probability distributions. By choosing 
A(ci,C2,(t) = (5(|g| — (?) we obtain the probability Pr{g) that two colliding particles have a relative speed 
|ci2| = 9- From here on we only quote results for two dimensions. Analytic calculations are based on the 
molecular chaos assumption (0) in combination with (0). Inspection of (11) shows that under this assumption 



the coUisional averages are independent of the x~f£^ctor. Straightforward algebra gives for the constrained 
(/—distribution, 

Pr{g) = {5{\c,2\-9))con 

Ve-^5'|l + la2(/-8g2 + 9)}. (21) 

Similarly we obtain the probability distribution for the center of mass velocity, G = i(ci + C2), 



PcM(G) = (5(|C|-G))con 

= 4Ge^2G^{l + a2(G4-G2)}. (22) 

It equals the unconstrained equilibrium distribution function apart from a small term of C(a2). Furthermore, 
the probability that the precollision speed |ci| of one of the colliding particles {i = 1, 2) has a value v is 

PH = (^(|ci|-i.))con 

= V2ve-'^'/mi+v^)Io {^v^)+v'h (^^)} , (23) 

whereas the unconstrained distribution is ~ t'exp(— u^). In evaluating this coUisional average we neglect 
the 02 contribution, and carry out the constrained <t integration. To calculate the remaining integral 
/dc2Ci2</5(c2), we change integration variables to g expressed in polar coordinates {g,(j)}, and use the rela- 
tion JJ^ d(?!)exp(— 2ci.gcos0) = TTlo(2cig). The subsequent g— integration follows from (6.618.4) in Rcf. p3^ . 
Using the asymptotic expressions /o(a;) ^ -^1(2;) ^ exp{x) / V^ttx for the modified Bessel functions of the zcroth 
and first order, we obtain the high energy behavior, P{v) ^ ^^f^pKv^ exp(— t;^). 

In a similar manner we obtain the following velocity moments and correlations, using the molecular chaos 
assumption, 

(5')coll = 3{1 + ifls}, (G2)eoU = Hi + 5«2}' 



4 "^z J 1 \^ /coil 2 



(.9*')coii = (g')con-(l-«'){l + ia2}, 



(c?)coii = (G2)co11 + 3(g')con = |{1 + ^02}, 
(ci ■ C2)coii — (G )coii —4(9 )con = ^jjl ^ 4^2}, 
(ci • c;)coii = (ci • C2)con + 5(1 - "^){1 + iaz}. (24) 

Here c* are the postcollision velocities, as defined in paper I. The sum of the third and fourth equality depends 
only on the center of mass velocity, i.e. (G^)coii- In the elastic limit a ^ 1, the average energy of a particle 
that is about to collide, (c^)coii = (5/4) (c^), is above the mean energy per particle, (c^), which equals unity. 

In the molecular chaos approximation an average like ((ci • C2)'"5")coii with {m,n} integers, is in gen- 
eral non- vanishing, except in the special case n = —1. Then ((ci • C2)™/(7)coii reduces to an unconstrained 
average, proportional ((ci • 02)"*), which vanishes for odd values of m. Additional information about the rel- 
ative orientation of the incoming velocities can be obtained from the distribution of the angle '012, defined by 
Ci • 02 = ciC2COS'0i2- A numerical calculation (again neglecting 02 corrections) gives (cos-0i2)coii — —0.233, 
which is close to the value —0.2, estimated from (ci • 02)0011 — (ci)coii(cosf/;i2)coii using the above results. 

A very sensitive probe for studying the violation of molecular chaos is the probability distribution P{b) of 
impact parameters, b = \g x a-\ = sm9, where 6 — cos^^(g-o-) is the angle of incidence. It is defined as the 
collisional average 

p.,^ /,,, i^^^n. /rf^/rfOi/d02(5(6-|gX^|)|g.^|e(-g-^)/(^)(ci,02,^) 

J d(T J dcij dc2 |g-cr|e(-g-(T)/(-^)(oi,02,cr) 

and P{b) can be easily computed in a molecular dynamics experiment. As long as molecular chaos holds, the 
distribution of b is independent of the functional form of / and we obtain straightforwardly 

Pib) = \^'-]i''-''[^<'<' (26) 

^ ' I otherwise, ^ ' 

which reduces in two dimensions to the uniform distribution, 

^ ' 1 otherwise. ^ ' 

In order to analyze molecular chaos breakdown in more detail, we have introduced a collection of moments 
M„„i and their dimensionless counterparts Bnm for n,m ~ {0,1,2...}, (see definition below), to analyze in 



detail the possible breakdown of the molecular chaos factorization (^. These moments Mnm{T) of the pair 
distribution at contact are defined as, 

M„„(T) = ^ydvidv2yd^e(-Vi2-^)/(2)(vi,V2,(T|T)<2|cos0r, (28) 

which are averages over the precollision hemisphere, where 6 — cos^^(g • a). Let A/^,„(T'e) denote the same 
quantity evaluated in Enskog's formulation of the molecular chaos approximation, and evaluated at the refer- 
ence temperature Te, i.e. evaluated with /(^^(vi, V2, o-|T) replaced by xe/(vi|Te)/(v2|Te), then the reduced 
moments are defined as 

Bn,n{T) = Af„™(r)/A/^„(rE), (29) 



where M^^{T-e) is evaluated in (A^). It is proportional to u^, where we = \/^2T^Jrn. We prefer to normahze 
the reduced moments by M^„(Te), because its analytic form is given explicitly. One could also normalize by 
M^^{T) = {T^/TY/'^M^^{T^). The disadvantage of M^^{T) is that the computation requires the simulated 
values of the kinetic temperature T. The coUisional averages (u"2l cos6'|'")coii expressed in terms of these new 
moments give, 

(Wi2|COS6'| ;coll - ^ /yx ■ (•aU) 

We first observe that the average collision frequency uj, defined in (||), is proportional to Afii(T), so that 

with We defined in (p9[). This implies that the reduced moments i?„„j(r) can all be expressed in coUisional 
averages, i.e. 

w W2"Vos^l"'^)con .„„. 

C^E W2 |COS0|™-1),E^;J 



The average (• • •)^q[j is defined through (30) with Mnm{T) replaced by Af„„j(TE), and calculated in (A. 2). It 



represents the coUisional average, evaluated with the Enskog factorization f^"^' = X'eI f o,nd also taken at the 
reference temperature Te. Note that the equality ( p2[ ) consists of two factors, u and (• • ■)coii, which are measured 
separately in event driven MD simulations. 

We also observe that the equality T{T) = rE(TE), explained above (p^, implies that 

Furthermore we have for the excess pressure, p°^(T) = p{T) — nT, 

„ (rj., pnT)/nT ^ oj {\vucos0\),,n 
^^^ ^ p'i{T^)/nT^ CE (bi2Cos0|)E^„ l^ ) 



(-) 



and for the dynamic pair correlation at contact x 

R lrr\ ^'"^ (|t>i2COs6l|-^)coll , - 

XE (|wi2Cos6i| i);^^n 

In the appendix we present a more complete set of relations for the Bnm ■ 

In the next section MD simulations will show that the predicted deviation from a Maxwellian (see (0)) in 
the (unconstrained) velocity distribution of a single particle can be observed for small inelasticities. However 
larger deviations are found between the observed constrained probability distributions and averages, and the 
corresponding kinetic theory predictions given by (|24|), based on molecular chaos. Consequently the small 
corrections resulting from 02 in (0) can be neglected in most cases. 



D. Radial distributions 

The static or unconstrained radial distribution function in the spatiaUy homogeneous IHS fluid is defined as, 

5(r)=y'^y'dcidc2/(2)(ci,C2,r^). (36) 

It may be averaged over all directions of r because of statistical isotropy. The unconstrained radial distribution 
function at contact is defined as the extrapolation, Y = g{r — > cr + 0). By splitting the it— integration into a 
precoUision (g • <t < 0) and postcoUision hemisphere (g • it > 0), we obtain Y as sum of two terms. 



Y 



i(r(-' + r(+)). (37) 



The definitions of F^^-*, y(+) follow from (B9) by adding respectively factors 0(— g • a) and 8(g • a) under the 
integral sign in (|3q). The dummy integration variables in y(+) represent the postcoUision velocities, (0^,02), 
corresponding to the precoUision ones, (01,02). 

On the other hand, we have the dynamic precoUision correlation x i defined in (|^), and a similar postcoUision 
one, x^^\ defined by replacing 0(— g • <t) in (ph by 8(g • a). They are related by continuity of flux. Because 
the incident flux of (C1C2)— pairs just before collision is equal to the scattered flux of (c|c2)— pairs just after 
collision, we have inside dynamic averages the equality, 

e(-5„)|5n|/^'nciC2,^)dcidc2d^ = e(g:)|g;|/(2)(cJc;,<T)dcJdc;da-, (38) 

where gn = g • o" = gcosO. The reflection law, g* = a|(7„|, for inelastic collisions in combination with the 
continuity of the flux and (H) yields at once. 



X 



(+) ^ n /ri.w(-) 



{l/a)x^-\ (39) 



In principle, equations ( |37| ) and (|20| ) provide two alternative routes to compute the precoUisional pair correlation 
at contact: the first one, denoted by Y^~' , can be implemented numerically as a static or unconstrained average, 
namely by extrapolation to contact of the pair correlation function for pairs aiming to collide. The second one, 
denoted by x^~\ can be computed as a dynamic coUisional average, calculated from /'^^^(ci, C2,(t) at contact. 
It is important to stress that the dynamic x is calculated as a time average over the subset of colliding 
pairs at contact, and the static Y^^' as a time average over all pairs, satisfying the relation, g • <t < and 
extrapolated to r -^ a + 0, i.e. calculated from /^^^ (ci , C2 , r) , where the limit is taken after all integrations have 
been performed. This may lead to different results, because the integrand contains the the function f^^^ which 
turns out to be singular near r = a and V12 small (see discussion in subsection III C). 

Physically, it is also clear why the averages in the NESS need not be the same. For instance, the relation 
may not hold for the limiting (r -^ a) values, y(") and y(+\ because non mean-field effects (in particular, the 
'rotation-induced' recoUisions discussed at the start of section III, or noise induced recollisions, see below) may 
result in differences between the two methods to evaluate x'^' and Y^~\ The reason is that the validity of (Bq), 
expressing flux continuity for the limiting values (r —> cr + 0), is questionable in the presence of the external 
random force. When the kicking frequency is much larger than the collision frequency (situation considered 
here), a pair of particles may indeed be put in contact under the action of the random force only. We will 
investigate possible numerical differences between x^~^ and y(^) in the next section on MD simulations. 

In paper I we have studied the long range spatial correlation functions Gabir) of the density field n(r, t) and 
the flow field u(r, i) in the NESS. These functions are closely related to (|3q), i.e. 



i j 

-Sir) + (gir) - 1) , (40) 



1 

n 



and, in the notation of paper I, 



-G||(r) + (d-l)G^(r) 

= ||<5(r) + «§g(,)(ci.C2)(r), (41) 

where (• • •) is an average over the A^-particle non-equihbrium steady state and the static velocity correlation, 
(ci • C2)(r), is defined as, 



(ci ■C2){r) = Jf dcidc2/(^)(ci,C2,rCT)(ci • C2)/.g(r). 



(42) 



The correlation functions Gabif) above are very long ranged, decaying like r^ '^ for large distances. In the first 
part of this subsection we have introduced the static correlations, Y ^ Y^^\ and the dynamic ones, x • In. 



(A. 8) and (|A.9|) we have done the same for the dynamic counter parts (ci • C2)dyn and (ci • C2)Jj ^ of the static 



. . . . ^_. - dyn 

correlation (ci • C2)(r — > a), introduced in (}41|). 

In the next section the short range behavior of these functions will be studied by MD simulations. 

III. SIMULATION RESULTS FOR THE NESS 

To investigate the short scale structure characterizing the NESS and the validity of molecular chaos we will 
present in this section MD simulation results, and compare these with our theoretical predictions whenever 
possible. The details of the simulations of the randomly driven inelastic hard disk system have been reported 
elsewhere pO| . We will work in the limit in which the kicking frequency of the external random force is much 
larger than the collision frequency. This is the limit in which the Fokker-Planck term in (Q) models the random 
energy input through the random kicks. The external random force will, in principle, have a quantitative 
influence on the short range structure of the fluid. There is only one important difference with respect to the 
simulations carried out in |1C| ] . There, the random rotation proposed in [|36| was implemented to avoid inelastic 
collapse at high inelasticities {a < 0.5). This procedure amounts to rotating the relative velocity g by a small 
random angle after each collision. Consider the completely inelastic situation a — for the sake of the argument. 
After each collision, the vector g lies exactly at the border of the precoUisional hemisphere (g-o- = 0), so that 
if the aforementioned random angle has equi-probable positive and negative values, the rotation procedure will 
lead to a recoUision with probability 1/2. This leads to a spurious increase of the number of collisions by a factor 
S^o 1/2" = 2 (the recoUision can itself induce a recoUision with probability 1/2 etc. . . so that the frequency of 
collision effectively doubles !). When a is small but non vanishing, this effect is still present but weaker. This 
is clearly an artificial violation of molecular chaos that has been discarded in the present work: for a < 0.5, 
we have also implemented the rotation procedure, but if a small rotation leads to a recoUision, a new angle is 
drawn until the pair separates. In this way, we reduce an important source of correlations (the effect is dramatic 
on all the low order moments Bnm, not only on the collision frequency; in particular, the moments with n < 1 
that correspond to coUisional averages of negative powers of g, are strongly biased toward bigger values if the 
"rotation-induced recoUisions" are present). After applying this new rule, we are then left only with correlations 
induced by the hard sphere dynamics plus the ones induced by the noise itself (see below). 

A. Cumulants 

First we focus on the single particle velocity distribution function / averaged over all particles, which deviates 
from a Maxwellian distribution due to the inelasticity of the collisions. In the previous section we presented 
predictions for these deviations, assuming molecular chaos. The resulting expression given by (0) for the fourth 
cumulant 02 of the distribution as well as the prediction for its overpopulated tail are in perfect agreement with 
3D Direct Simulation Monte Carlo (DSMC) results over the whole region of inelasticities 134]. As DSMC itself 
invokes molecular chaos, this observation merely justifies the approximations made in the analytic calculation. 
Information about the validity of molecular chaos can only be obtained from a comparison with molecular 
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dynamics simulations, and Fig. g shows this comparison for the fourth cumulant 02 as a function of the 
coefficient of restitution a in 2D. The simulation results are in agreement with (Q) for small inelasticity, but 
start to deviate significantly from the theoretical prediction at a = 0.6. These deviations, together with 
the perfect agreement between the theoretical prediction and DSMC results, provide direct evidence for the 
breakdown of molecular chaos for a < 0.6. The theoretical result (0) is independent of the density. As a^ 
represents only a small correction in (y) , one needs a large number of collisions and a large number of particles 
to reach sufficient statistical accuracy. So, high densities, for which one can use linked lists [Q, are well suited. 
The data in Fig. g are typically obtained at high packing fractions (0 = 0.63 at a = 0.92 and 0.7, and = 0.55 
at a = 0.6) and N — 10201 particles. At low densities and weaker inelasticities (a > 0.6) we are unable to 
collect enough statistics to measure the small correction, represented by a-i. Simulations at higher inelasticities 
did not show any density dependence of a^ in the range 0.2 ^ 4> ^ 0.6, suggesting that the cumulant expression 
(ra), obtained from the Enskog-Boltzmann equation also applies at liquid densities. Similar results as those 
displayed in Fig. have been observed for the 3-dimcnsional version of the present model pT[ |. 

B. Radial distribution function 

Next, we present results for the static or unconstrained radial distribution function g{r), and in particular its 
extrapolated values at contact Y, Y^~\ Y^~^\ Here g{r) is essentially the density-density correlation function, 
whose long range behavior was studied in Ref. |lQ| . 

Figure g shows the measured values of g{r) for short distances and packing fraction cj) = 0.2, at different 
inelasticities. At small inelasticities (a ~ 0.9), g{r) resembles the radial distribution function for elastic hard 
disks (EHD). At higher inelasticities, deviations start to appear: the first and second maximum in the measured 
g{r) are enhanced with respect to their EHD values at the same density. Moreover, the functional shape also 
deviates from the corresponding pair distribution of EHD at an appropriately chosen higher density; e.g. if this 
density is chosen such that the value of the second maximum of the pair distribution of EHD coincides with 
the simulation result for IHD, the observed value at contact would still be underestimated by the EHD pair 
distribution. 

It seems worthwhile to compare these results with existing experiments on granular fluids in which the pair 
distribution g{r) has been measured. In the experiment of Ref. |3| on a vertically vibrated thin granular layer, 
g{r) has been measured a.t <j> — 0.46. In the fiuidized ('gas-like') phase, it follows the equilibrium result for 
elastic hard disks almost identically. This result may be compared to our simulations for a randomly driven 
fluid of inelastic disks at a = 0.9, corresponding to the value for stainless steel balls used in the experiment. 
It would be of interest to measure experimentally how g{r) in the fluidized phase depends on the inelasticity, 
and see if a behavior similar to that of Fig. pi is observed. It is also interesting to note that the pair correlation 
function g(r) in a non-Brownian suspension of spherical particles, fluidized between two vertical parallel plates, 
shows an enhanced value at contact as well p^. 

In Fig. Ha we show the value at contact, Y, obtained by extrapolation from g{r) at = 0.05, together with the 
extrapolated values for approaching and receding pairs, y(~) and y(+' respectively. For a ^ 0.8 no significant 
deviations are found from the Verlet-Levesque value xe = 1.084 for elastic hard disks at the same density. More 
surprising is the value oi g{r = a) at large inelasticities, reaching a value around 40 for a — > 0. This property, 
combined with the observation that the first and second maximum in g{r) are shifted to smaller r-values, and 
are larger (up to 20% at small a) than the corresponding hard disks values, may be interpreted as a tendency 
to cluster, i.e. to stay in continuously rearranging configurations with large density inhomogeneities. We return 
to this point in subsection HI G. 

Figure also shows the dynamic correlation, x = Xe-Bqo, measured as a collisional average. Fig. 
compares the static ratio y(~)/y(+) with the dynamic one, X^~Vx*'^'' = 0^1 ^ind also shows the ratio l''~Vx^ 
The plots clearly show that the dynamic and static correlations, x^^^ and Y^^\ are different. For a ^ 0.5 the 
differences are large, and for a > 0.6 both functions are about equal. For the case of a freely evolving IHS 
fluid, Soto and Mareschal [|8| have recently observed a similar behavior, and explained it in terms of the effect 
of the increase of grazing collisions on the effective x^~'*- In the randomly driven IHD fluid the same effects 
are present. All correlations, x''~'*i l"*-^-* are large, especially at small a. This is caused by the divergence of 
/'-^-'(ci, Ci, cr) at small g and small cos 6*, which corresponds to grazing collisions and will be further discussed 
in the next subsection. As a result of noise-induced recoUisions, collisions with small g and small cos9 are 
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oversampled; consequently, a dynamical average involving negative powers oi gcosd such as x is expected to 
be larger than its static counterpart Y^~\ This feature can be observed in Fig. |jb. Moreover, in the absence of 
recoUisions, we would expect y(~) = aY'--^^ as a result of plain hard sphere dynamics. However, in the heated 
system, the flux continuity, as expressed in Eq. (p8|), is no longer satisfied for the extrapolated Y's; some pairs 
are put in pre-coUision configuration under the action of the random force, which leads to Y^~' > aY^^' . The 
breakdown of Y^^' = aY^^' signals the inelasticity beyond which noise-induced correlations become relevant. 
It is furthermore possible to consider situations where the recoUisions dominate the dynamics, e.g. at small a 
by allowing rotation-induced recoUisions. In this extreme limit, we expect Y^^' ~ Y^'^' , the pairs being put in 
pre- or post-collision configuration essentially at random. In the same limit, the population of colliding pairs 
with small g and small cos 6 is enhanced, leading to a more pronounced discrepancy between dynamic and static 
averages (i.e. a much smaller ratio Y'-~^ /x^~^ than observed in Fig. 0b). 

C. Equation of state. Molecular chaos breakdown 

To what extent does the extrapolated static radial distribution function, Y = g(r —f a), describe the nontrivial 
dependence in the NESS of collision frequency in (^, collisional damping in (|lj) and pressure in ( [l2|) on the 
inelasticity? If molecular chaos holds, the latter quantities depend, according to Eqs. ffl), (O), (H) and (llq), on 
the prccollisional pair function at contact, x ^ where the particles are aiming to collide. This function differs 
from the extrapolated static Y'^ ' at high inelasticities (see Fig. H). Consider first the collision frequency in the 
molecular chaos approximation, o^mc ~ X'^o{T) above (|9|), with x ~ X the dynamic correlation, i.e. 

r>00\ TfT^ I4dj 



^E Xe ve V Je 

where we have used Eqs.(Bf) and (pq). This is an extremely poor approximation, as can be seen from Fig. m, which 
shows that the measured value lJ/uj-e approaches 5.6 as a —> 0, whereas Bqq is essentially divergent. Next we 
replace x^~^ in ( ^3| ) by its static counterpart Y^~'\ shown in Fig. ^. This yields iVstatiT) /uje = Y^~^vo/x'"~^''^e- 
Its limiting value for a — > is about a factor 3 too large when compared to lo/uje- We conclude that all mean 
field approximations for the collision frequency, including the Enskog approximation u!e{T)/uj-e = yT/T^, 
break down for a < 0.6. 

Fig. H shows the pressure of the IHD fluid, compared with the molecular chaos prediction given by (13), taking 
for X either the Enskog approximation xe in (|3)j or the simulated Y^~' , or x ■ The Enskog approximation, 
accounting for the short range geometric exclusion effects in the precoUision state, gives a reasonable description 
oi p{T) for all a, while both the static F(~) and the dynamic x^~^ give an extremely poor description except 
for a > 0.8. 

Consistent with this conclusion is the good estimate for the temperature Te in the NESS, obtained by 
balancing the energy dissipation rate rE(rE) in (Hq) with the energy input from the random force, as shown 
in Fig. 0. Moreover the colUsional energy loss, TIT)/Te(T) = {T^/Tf^'^ in (^, is in agreement with MD 
simulations over the whole a- interval within 30%. All other mean field approximations with t^E replaced by 
Wmc(T) or ujsta.t{T) give very poor results for T[T) = mQ. 

How can these paradoxical results be reconciled? Let us compare the individual definitions of x t^tP ^-nd 
r, which all contain factors |g • cr|"/^^^(ci, C2, <t) with n = 0,1,2,3. To find a possible explanation of these 
paradoxical results, we test the following scenario: the molecular chaos assumption (|gj only breaks down at very 
small relative velocities g, and more precisely, at very small gn ^ S ■ o' = gcosO, which is the component of g, 
parallel to the line of centers of the colliding particles (physical arguments for this scenario will be offered in 
subsection IIIF| where we discuss the noise induced recoUisions). On the basis of this scenario the singularity 



in /(^' at small g makes the dynamic correlation x — ^ooXe ( shown in Fig. p) very much larger than 
Xe, essentially divergent as a — > 0. In calculating the collisional frequency from (ra) the extra factor g„ in the 
integrand makes the small o„— singularity integrable, giving a finite correction to the Enskog collision frequency, 
also for a ^ (see Fig. |^). The contributions of the small g„— singularity in p{T) and r(r) are essentially 
suppressed by extra factors of |g • (t|". 

This possibility has been analyzed systematically by measuring the behavior of the moments Bn„i{T), which 
are useful tools to investigate the breakdown of the molecular chaos assumption. We have made n and/or m 
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small in order to analyze the nature of the singularities in f^'^^ near small relative velocities and near grazing 
collisions, as displayed in Fig. |7|. All deviations of these quantities from unity give a quantitative measure for the 
violation of the molecular chaos assumption. In the elastic limit we have carefully checked for a large number of 
cases that the reduced moments Bnm tend to 1. Fig. g shows the values of different moments Bnn(T), and one 
can clearly see how the deviations from the elastic limit rapidly decrease as n increases to n = 3, after which 
they start to increase slowly. For larger n-values, the moments are reasonably close to unity, but statistical 
inaccuracy precludes any definite conclusion about the large n behaviour. 

Further evidence for the above scenario is shown in Fig. |7|, where we display two sequences of moments Bnm- 
To draw some further conclusions from Figs. ^ and |^ we note that the integrands in Bnm Bon and Bno, as 
defined in (psj), contain apart from /^^^ respectively the factors g"\ cos0|", | cos9\'^,g^. The reduced moments 
Bqi and Biq contain again very large contributions from the divergence of f^'^' near vanishing g„ = gcosO. 
Fig. suggests that the presence of equal powers of g and cos 9 in i?„„ simultaneously suppresses the large 
contributions from the singularities at g = and cos 9 = Q. 

We conclude that the numerical results, displayed in Figs. ^ and ^ give support to the previous scenario, 
showing that molecular chaos breaks down only in a very small portion of the phase space, around gn = 
g • CT = gcos9 = 0. The size of this 'pocket' in phase space increases as a decreases. Therefore, only those 
coUisional quantities that contain low powers of g and cos 9 (such as x*-""* and uj) will be very sensitive to this 
breakdown as the inelasticity increases, while physical quantities involving higher powers of g and cos0, such as 
the temperature, pressure or energy dissipation will be well approximated by their molecular chaos counterparts. 

D. Velocity correlations at contact 

In the previous subsection we have considered the pair distribution function /"^(ci, C2, cr) in the precoUision 
state, and have examined how molecular chaos is broken down, and which physical quantities are most sensitive 
to it. Now we will analyze the effect of the breakdown of molecular chaos on coUisional statistics. 

We show in Fig. pi different velocity coUisional averages at cj) = 0.05. In the simulations, these quantities are 
obtained by averaging over successive collision events in the steady state. We first observe that the simulation 
results in Fig. ^ approach for a — > 1 the analytic results for elastic spheres, calculated in (p4|). At small 
inelasticities, the simulation data follow the trends of the theoretical prediction with systematic deviations 
depending on the quantity considered. For instance, the behavior of the center of mass motion (G'^)coU is 
close to the analytical prediction of ( p4[ ) in the whole range of a values. This indicates that the center of 
mass velocity G is not correlated with the relative velocity g. Consequently, /(^)(ci,C2, cr) in the coUisional 
average (O) factorizes, and we can expect the contributions in numerator and denominator in (0) coming from 
G— integrations to cancel. Consistent with this behavior, we observe that the two curves in Fig. M , (c^)coii 



(labeled by circles) and (ci • C2)coii (labeled by squares), are symmetric around 1/2. In Eqs. (A. 4) and 
the appendix these quantities have been expressed in reduced moments. 




(Cl • C2> 



1 3631 

2 4 6n 

1 3&31 

2 4 6ii' 



(44) 



where 



.(r) fT, 



M^„,{T) V T 



tE 



n/2 

Bnm- (45) 



The reduced moments have been measured independently (see Figs. and R), and used to calculate the 
expressions ( [44| ) and (^). The results have been plotted in Fig. ^ as dashed and dashed-dotted lines, which 
agree very well with the direct measurements of these quantities as coUisional averages, shown in Fig. 
respectively as squares and circles. In deriving (|4J) and ( [45|) we have again used that the velocity variables G 
and g are statistically uncorrelated. The present results strongly support this assumption. 

The correlation (ci -£2)0011 — (cos V'i2)coii, also plotted in Fig. H, cannot be expressed in B-moments. However, 



the approximate relation already employed to show that (cos?/'i2)coii — (ci • C2)coii/(c^)coii in section II C, holds 
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for the simulation data over the whole range of inelasticities. As the system becomes more inelastic, the typical 
"temperature" of colliding particles (defined as the coUisional average (c^)coii) decreases and even becomes 
lower that the unconstrained average (c^) that defines the temperature. On the other hand, as already noted 
below (gj), (c^)coU = 5/4 > 1 in the elastic limit. This decrease of (c^)coii is directly related to the increase of 
the small gi-portion of phase space where molecular chaos is violated. At small a, most of the collisions occur 
between particles with small and even vanishing relative velocities. An extreme example is the inelastic collapse, 
mentioned in the introduction. 

The correlation function (ci • 02/5)0011 for the freely evolving IHD fluid has been simulated by Soto and 
Mareschal, and was shown to be small, but non- vanishing pq |. 

We have also investigated r — w— correlations by measuring the expectation value of (ci ■C2)(r) for two particles 
separated by a distance r, as defined in (p^). The results are shown in Fig. S and Fig. M The plot shows an 
intermediate range of r values with an exponentially decaying correlation. It is again of interest to compare the 
extrapolation of the static correla tion (ci • C 2)(r — > a) with its dynamic counter part, (ci • C2)dyn, calculated at 



collision. The results, derived in (|A.§|) and (A. 9) of the appendix, read for hard disks. 



/ \ 1/^1 ^2o\ 1-a ^22 ,,„. 



The first term on the RHS represents the precoUision part. 



(ci • C2)U = Ul - ^ ) ■ (47) 



uoo 

Figure |0| compares the extrapolation (ci • 02) (r -^ <j) (circles) of the static correlation with its dynamic analogs 
(EgI) and (p7|). The numerical data for both correlations agree well for a ^ 0.8, but for a < 0.5 the dynamic 
correlatioii(solid line) is substantially larger than the static one. This is consistent with the difference between 
X^"-* and y(^) observed in Fig. ^. For comparison the dynamic precoUision correlation (dashed line) is also 
shown . It should be noted that the divergence of J*^^^ at small g and small cos 6 implies in particular that 
Bqo '3> B20 > -B22, so that ( [l6|) predicts that the dynamic correlation at contact (ci • C2)dyn should increase 



at a — > and saturate close to 1/2. By the same arguments its precoUision part in (47) approaches the same 
limit. This can be observed in Fig. |l^. 

The velocity correlation (ci • 02)0011 in (M) involve the reduced moments 631 and 611. Consistent with the 
scenario, developed in subsection III C, the divergence of /(^^(ci, 02, u) near g — and cos9 = is largely 
suppressed in these higher moments, which remain finite for a — > 0, where 611 ~ 4631. Consequently (ci • 02)0011 
does not approach the value 1/2 as a ^ 0, but a value close to 0.3, as can be deduced from Fig. g. 

E. Grazing collisions 

The data in Fig. ^ for (oi • 02)0011, (cosf/;i2)ooii, (Vl — b'^)coU and (6)ooii clearly illustrate that the violation of 
molecular chaos strongly increases with increasing inelasticity. Consider first the average, 

(&)ooii= / dbbP{b). (48) 

Jo 

This average remains at a plateau value 1/2 for a ^ 0.5 , which is determined by the uniform distribution 
P{b) corresponding to molecular chaos in two dimensions. Recall that the value 1/2 holds regardless of the 
functional form of the velocity distribution function /. It is thus a good probe for molecular chaos breakdown. 
Moreover, from its trend we can also estimate the way in which such a breakdown takes place. Specifically, 
as the inelasticity increases the average value increases by about 50%, which indicates a strong bias toward 
grazing collisions. To illustrate this, we model the normalized distribution of impact parameters as a uniform 
background and a 'half delta-peak a,t b = 1, i.e. P{b) = 1 — p + 2p6{l ~ b), where p is the fraction of grazing 
collisions. This yields the average (6)coii — 5(1 +p), which implies, according to Fig. g, that at a = 0.0, 0.1 and 
0.3 respectively a fraction of 50%, 35% and 5% is grazing at = 0.05. This qualitative picture is supported in 
a more quantitative manner in Fig. nil which shows the measured P{b), which is strongly peaked near grazing 
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collisions (6=1). At small inelasticity all impact parameters are equally probable as expected on the basis of 
molecular chaos, and consistent with Fig. ]q. Only for a ^ 0.5 deviations become significant: upon decreasing 
the coefficient of restitution, collisions with a larger impact parameter occur more frequently, implying an 
increase of the frequency of grazing collisions. The behavior of P{b) is then fully consistent with the divergence 
of /^^^ at small cos6', discussed in subsection III C. 

To avoid inelastic collapse for a < 0.5 the postcollision velocities of colliding pairs are rotated over a small 
random angle as described in Refs. |pq,Qffl, with the important restriction mentioned at the beginning of section 
III. Alternative algorithms to avoid inelastic collapse are described in Ref. |g^. For a > 0.5 no such rotation 
was applied. To check if the deviations of the impact parameter for a < 0.5 are due to this applied rotation, 
we have also performed simulations where even for a > 0.5 a random rotation was applied. Regardless of the 



applied random rotation, we found (6)con close to ^ for a > 0.5. Both Figs. H and 11 show that for a ^ 0.5 
molecular chaos is strongly violated, and that the violation is weaker in the small inelasticity regime. The 
average (Vl — 6^)coU supports the same conclusions. 

The data for (ci • 02)0011 and (cos'(/'i2)coii in Fig.0 are consistent with the predominance of grazing collisions 
at large inelasticities. They show the average relative angle between the velocities of the incoming particles, 
which has a strong a— dependence and no plateau value near the elastic limit. Near a = 1 the particles are on 
average on approaching trajectories with (cos -012)0011 — —0.25 and (i/'i2)ooii — 105°. As a decreases, (cos -012)0011 
increases linearly to a value 0.50, while (^12)0011 approaches 60°, at a = 0. This corresponds to collisions of 
more or less parallel-moving pairs of particles, where faster particles overtake slower ones. 

Figs. O and O show the distribution of relative orientations of incoming velocities. The distribution of 
angles between the incoming particles (V'12) shows moderate deviations from what is expected for an elastic 
system in the range 0.5 '^ a < 1. As an analytic expression for elastic disks is not available, deviations are 
compared with the simulation results for elastic hard disks (in the absence of a random external force). At 
a = 0.5 the frequency of collisions of parallel-moving particles is strongly increased, a trend which is enhanced 
upon increasing the inelasticity. Finally, the probability distribution P(ci • C2) is shown in Fig. O. When 
the inelasticity increases, this distribution becomes more peaked around the origin, as the colliding particles 
on average move more slowly relative to each other. In the mean time, the typical angle ^12 decreases, which 
causes this peak to shift to positive values. 

F. Particle- and noise-induced recollisions 

The mechanism for the breakdown of molecular chaos in classical fluids with conservative interactions are 
sequences of correlated ring collisions, as discussed in the introduction. The most simple ring collisions are the 
recollisions (1-2) (1-3) (1-2) and cyclic collisions (1-2) (2-3) (3-1) or permutations thereof pO[ |. 

There is strong evidence that the effects of ring collisions are considerably enhanced in fluids with dissipative 
interactions, such as granular flows, where relative kinetic energy is lost in binary collisions. As a result the 
postcollision velocities {v^, V2} will be on average more parallel than the precoUision ones {vi, V2} pi], i.e. the 
trajectories are less diverging than in the elastic case, and there is a much larger {rs, V3} phase space, in which 
particle 3 will knock, say, particle 1 back to recoUide with particle 2. 

This increase of phase space is confirmed by gathering recollision statistics. We have counted the fraction of 
recollisions as a function of a, as shown in Table Q. The column labeled 7?,i (recollisions between two partners 
mediated by a third particle) shows that at a packing fraction (p — 0.2 in the elastic case {a — 1) only a fraction 
of 6.7% of all collisions is a recollision. This frequency gradually increases to about 15% at a = 0.4. 

In the randomly driven IHS fluid there is the additional effect of noise induced recollisions that do not require 
the intervention of other particles. This type of recollision (denoted TZq) occurs with high probability when the 
relative velocity after collision is so small that it can be simply reversed by a random kick. At a = 0.6 the 
frequency of noise-induced recollisions is about 4%, and it increases to 52% at a = (see column TZq in table S). 
The effect is of importance at all densities, because it does not require the mediation of a third particle. Indeed, 
at a low packing fraction of 1% and in the completely inelastic case a = 0, the frequency of 7?.o-like events is still 
34%, while 7?.i-like events have dropped to 5%. Moreover, we have verified that inclusion of rotation-induced 
recollisions modifies most of the collisional quantities we have analyzed, increasing their deviations with respect 
to the molecular chaos prediction. 
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At present, more quantitative theories or estimates of the effect of both types of recolHsions and other ring 
colhsions on the short range behavior of the pair distribution function /"'(zi, X2) are lacking. A natural way 
to incorporate the noise-induced recollisions into a kinetic theory description would be to include them into an 
effective two-particle scattering operator, which transforms an asymptotic precoUision state of two independent 
particles into an asymptotic postcoUision state, without involving intermediate two-particle scattering states, 
as in the present case. This may lead to an instantaneous Boltzmann collision term (without memory effects), 
provided the mean free time and the time between random kicks are very well separated (dilute gases). Such a 
description would suppress the recoUions of type TZq j and make the violation of molecular chaos less severe, say 
comparable to the freely evolving IHD fluid. 

G. Cold dense inhomogeneities 

In Ref. pc| ] we have shown by analyzing the Fourier modes of the granular hydrodynamic equations, which are 
valid for small inelasticities (say a > 0.7), that the NESS in a randomly driven IHS fluid is linearly stable against 
spatial inhomogeneities. Consequently, when observed over sufficiently long times, the NESS should be spatially 
homogeneous. However, it was also shown that the NESS exhibits strong fluctuations, resulting in long range 
spatial correlations in density, flow field and granular temperature. The observation of density inhomogeneities 
for large inelasticities has already been reported by Peng and Ohta [||. These density inhomogeneities, as shown 
by the snapshot of the density in Fig. nj, are not quasi-static, as in the freely evolving case |38,pq,p9,pJ| , but 
seem to behave as dynamic assemblies of particles that dissolve and re-assemble again. Also for a uniform shear 
flow, dynamical density inhomogeneities have been reported Effl. The existence of density inhomogeneities was 
already suggested by the static pair distribution functions g{r), which showed an enhancement of the first few 
maxima as compared to their elastic values (see Fig. H). 

In Fig. M we show that the mean energy (cf)coii, of particles aiming to collide, is above the mean, (c^) — 1, 
for small inelasticity. It decreases from its elastic value 5(c^)/4 with decreasing a, then crosses the mean value 
value (c^) = 1 at a ~ 0.2, and further decreases to approximately 0.7(c^) at a = 0. 

It is interesting to observe that in the strong dissipation range, the mean kinetic energy or granular temper- 
ature of particles that are about to collide is lower than the average temperature. We combine this observation 
with Figs. 3a,b of Peng and Ohta [BJ, which show that essentially all collisions occur inside "cold" regions of 



high densities. This last observation applies even more so to undriven IHS fluids |38 0|. We expect that, also 
in the randomly driven IHS fluid, the majority of collisions takes place inside cold high density regions. 

If the predominance of cold particles in strongly inelastic collisions, (c^)coii < (c^) is indeed a signal for 
the appearance of density inhomogeneities, then Fig. g suggests that at a packing fraction (j) = 0.05 density 
inhomogeneities may occur for a < 0.2. This is indeed confirmed by the snapshots in Figs. 111. In Fig. W^ 
we illustrate the existence of cold inhomogeneous dense regions for a — 0.2 and (p = 0.2. The particles with a 
less (more) than median kinetic energy are shown on the left (right). The formation of inhomogeneities is more 
clear for the colder particles. The temporal evolution of these regions show that they dissolve after some time, 
while new inhomogeneous regions appear. The formation of "living" inhomogeneous regions can be understood 
using the hydrodynamic picture put forward in llG], where it was shown that the structure factor behaves as 
S'(k) ~ k~^, implying density correlations increasing with distance as ln(r) in two dimensions. These long range 
spatial correlations induce a slowing down of the dynamics, as in critical phenomena. This, in turn, implies the 
slow decay of density perturbations, that could lead to visible density inhomogeneities as the kicking frequency 
is reduced (in this respect, see Refs. M). We can also expect that upon decreasing the forcing frequency, the 
dynamics should be closer to its "free cooling" counterpart so that well defined clusters are then likely to appear. 

More details about the predominance of cold particles, among those involved in collisions, can be seen in Fig. 
n6^, which shows the constrained probability distribution P(c), defined in subsection II C and obtained from 
MD-simulations at different inelasticities. For a ^ 0.5 the distribution has significantly shifted to smaller impact 
velocities. For the completely inelastic case, collision events involving "immobile" particles are more than twice 
as frequent as for the elastic case. The second moment of the distribution displayed in Fig. |l^ decreases when 
increasing the inelasticity. In fact, all functional forms with simulation data at different a can essentially be 
collapsed onto a single universal curve (the clastic one) by plotting ■y/T(a)P(c|a) as a function of c/ yjT{a), 
where T{a) = (c^)coii is the mean temperature of a particle at collision. The collapse plot is shown in Fig. [iqb. 
This data collapse confirms the concept of cold dense regions dominating the energy dissipation. This could 
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point to a possibly relevant two fluid picture of a "hot" dilute background gas coexisting with continuously 
rearranging configurations of "cold" dense regions. 

IV. CONCLUSION 

We have performed extensive MD simulations to study the kinetic properties and short range correlations 
in the non-equilibrium steady state of a randomly driven fluid of inelastic hard disks, as a model for fluidized 
granular material. The MD results have been compared with kinetic theory predictions derived from the 
Enskog-Boltzmann equation, properly modified with a Fokker-Planck diffusion term ^q(9/9v)^ to account for 
the applied random driving force P]. 

It appears that the kinetic theory predictions, based on molecular chaos, are essentially in agreement with 
the MD results for small inelasticities {a > 0.5) at ^ = 0.05. For larger inelasticities the deviations from the 
molecular chaos predictions start to become manifest: the radial distribution function at contact differs strongly 
from its local equilibrium form; there is a predominance of grazing collisions. When increasing (jj, the effects of 
the inelastic collisions become relevant at smaller inelasticities; e.g. a.t cj) = 0.2 and (j> = 0.5, we observe already 
significant deviations for a < 0.7. 

To avoid inelastic collapse of the system at low a, we have implemented a modified rotation procedure (see 
the beginning of section III). In its original version, this procedure induces dramatic violations of molecular 
chaos. It could then be argued that the important deviations of low order reduced moments Bnm are also 
spurious consequences of the above rotation procedure. However, we checked that circumventing the collapse 
by applying elastic collisions when the relative velocity of a pair is below a certain cutoff p7| ] , also induces very 
important violations of molecular chaos (quantified by -Boo for instance), unless the cutoff is chosen unphysically 
high. 

Sequences of ring collision processes, which lead to the breakdown of molecular chaos in classical fluids with 
conservative interactions, are strongly enhanced in fluids with dissipative interactions, like rapid granular flows. 
We have analyzed how molecular chaos is broken, i.e. essentially only through pairs of colliding particles at very 
small relative velocities. This means that molecular chaos is violated only in a small portion of phase space, 
implying that only certain physical properties will be sensitive to this violation. This explains why quantities 
like the collision frequency, or the pair distribution function at contact are very sensitive to the inelasticity 
parameters, while others like the pressure or the energy dissipation rate are well approximated by their Enskog 
prediction. Disentangling the effects of hard disk and noise-induced correlations remains an interesting point to 
explore. The studies performed in a freely evolving IHS fluid also shows the predominance of grazing collisions 
at long times. The fact that we have observed an analogous behavior for this homogeneous steady state indicates 
that the mechanism of breakdown of molecular chaos in granular fluids through grazing collisions is generic for 
this type of fluids. 

The extra feature of noise-induced recoUisions, which do not require mediation of a third particle, will further 
enhance the violation of the molecular chaos assumption. A natural way to develop a kinetic theory for randomly 
driven fluids, thereby presumably restoring the validity of the molecular chaos assumption in the dilute gas case, 
could be to include the noise-induced recoUisions in an effective two-particle scattering operator. It would be 
of interest to study its properties, either analytically or by simulating a two-particle inelastic collision in the 
presence of external noise. An additional theoretical complication here is the validity of the Boltzmann equation 
(p with Fokker-Planck diffusion term due to the fact that there are two limits involved when dealing with hard 
spheres in combination with external white noise. The actual properties of the effective collision operator depend 
on the order in which both limits are taken. In the simulations, one always takes the hard sphere limit first, 
while the white noise is approximated by discrete kicks which are applied to the particles at discrete times. 

In Ref. p0| we have calculated the equal-time spatial correlations of the fiuctuations in the hydrodynamic 
densities in the NESS. Here we have focused on the dynamic properties of these enhanced fiuctuations, in 
particular of the dynamic inhomogeneities observed in the density field. The coUisional velocity moments, 
introduced in section y and measured in MD simulations, reveal that the dense regions consist mostly of 
particles colder than average. This is clearly shown in the velocity distribution P(v\a) of particles which are 
about to collide. 

The MD simulations have been carried out in the limit in which the time interval between the external random 
kicks is much shorter than the mean free time between collisions. In this limit, regions with density larger than 
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average are not seen to survive for a long time. Rather, they form, dissolve and reappear elsewhere. The spatial 
correlations analyzed in |10(] show long-range correlations, which imply also a slowdown in the temporal decay of 
density perturbations. Therefore, we expect than the decrease of the kicking frequency will be accompanied by 
the appearance of apparent clusters. This fact, together with the shape modification of the velocity distribution 
P{v\a) (see Fig. nw suggests the picture of a two fluid model, in which a "hotter" more dilute background 
gas coexists with continuously rearranging configurations of "cold" dense clusters. This point remains open for 
subsequent investigation; for example, it would be interesting to analyze separately the coUisional statistics in 
the dense and dilute regions to assess the role of density fluctuations. 

APPENDIX: REDUCED MOMENTS Bnm{T) 

In the body of the paper we have considered the coUisional averages (^"l cos6'|™)coii and the moments Mn,m{T) 
and Bn,m{T)- We first list the Enskog values of these quantities, which have have been calculated from its 
definitions, given below (psf), i.e. 

AfE lTu\ - 7;"v,.^"/^ r((d+«)/2)r((m+l)/2) ..,x 

/«"lrn«fll™\E _ on/2 r((d+«+l)/2)r((m+2)/2) ..„n 

\9 |C0SW| ^eoll - ^ r((rf+m+l)/2) ' ^^■'^1 

Many physical quantities of interest can be expressed in terms of reduced moments Bnm , as already illustrated 
in subsection II C for x'^^cjjp and T. Analogous relations hold for the velocity moments (5")coii, which are 
proportional to M„+i^i. This yields 

^ (1^12)0011 _ w (.g")coii / r ^ "^^ 



B„ 1 = — X p-^^ = — X ppp — , (A.3) 



where the denominator has been calculated in (A.2) 



Velocity correlations between nearby particles can also be expressed in the reduced moments Bnm{T). First 
consider the constrained averages (ci • C2)coib defined in (|24|). They contain {G'^)co\l^ which equals d/4 from the 
MD simulations, in agreement with (24) (see Fig. |8| of section III). The center of mass velocity G is consequently 



uncorrelated with the relative velocity, and independent of the inelasticity. Substitution of (G^)coii = d/A: in 
HI) yields, 

(ci • C2)con = 4 - 4(.9 >con = 4 - ^ [y ) ;^ (A.4) 

/ 2\ ci ,1/ 2\ d d+1 /Te\ B31 

(ci)con -4 + 4(5 )con = 4 + ^ [-y) :^ • (A.5) 



Similarly we find. 



(ci-C2/.9)eon = (GVff),on-i(-9) 



/coll 



d Tid/2) [T ^B^(^(T_^\B^ ^^^^^ 



4v^r((d+l)/2)VTE Sii i \T J B, 
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(ci • C2/|gcos0|),,i, = (GVl.gcos0|)^^„ - \ (.g/| cos0|)^^i, 
d pn fr' Boa f, /TeX -B21 




and 



(A.7) 

DQl J 

Note that the last two averages are vanishing in the elastic case. 

In the body of the paper we have considered the extrapolation of the static correlation (ci • C2)(r -^ a). 
Here we calculate its dynamic analog (ci • C2)dyn, obtained by interchanging limits and replacing f^'^'{ci,C2,r) 
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under the integral sign in Eq.(|4l|) by its value at contact, /'■^^(ci, C2, o"). We proceed in the same fashion as in 



Eqs.(|36|)-(39), and split the numerator in ( [42[) in a pre- and postcollision part, as done in subsection II D. One 



finds after a lengthy calculation. 

Here the first term on the RHS is its precoUision part, i.e. 

(Cl •C2)|iyn = (Ci -0215 cos 6*1 "^>coll / (Iff COS 6*1 "^)con 

In section III these quantities are compared with MD simulations. 
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7^o 


7^l 





52 % 


18% 


0.4 


14% 


15% 


0.6 


4% 


15% 


0.95 


0.15 % 


7% 


1.0 


0% 


6.7% 



TABLE I. Frequency of recoUision events as a function of the inelasticity (see text for the definition of TZo and TZi) 
The packing fraction is = 0.2 and the system contains N = 5000 disks. 
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FIGURE CAPTIONS 

1. Kinetic temperature T/T^ and collision frequency uj/lu-e where Te and cje are defined in Eqs. ( |l7|) and 
(p^), for a packing fractions cf) — 0.05 and = 0.2. 

2. Fourth cumulant as a function of the coefficient of restitution. Comparison is made between the two- 
dimensional version of M) and MD results (circles) obtained for a system of 10201 inelastic disks, measured 
at several densities (see section III A). 

3. Pair distribution functions g{r) versus distance between the particles at a packing fraction (p = 0.2. The 
arrow indicates the value at contact for an elastic hard disc (EHD) fluid (from Verlet and Levesque ^^). 

4. a) Static or unconstrained pair distribution functions at contact Y,Y'-^\Y^^'' , as extrapolated from the 
corresponding pair distribution function g(r — *■ cr) at ^ = 0.05, compared with the dynamic correlation 
X^^^ at contact. The straight line corresponds to the EHD prediction (xe — 1.084). b) Ratio of dynamic 
to static correlation x''~V^' ') to be compared with 1, and the static ratio Y'--^'/Y^^\^ to be compared 
with the dynamic ratio X Vx = Oi- 

5. Pressure versus coefficient of restitution at a packing fraction cp = 0.05. The simulations results (direct 
or through -B22) are compared with molecular chaos prediction (|l3|), where x is either the static y(^) in 
(|37|), or the dynamic x in (pD? or the Enskog approximation xe in (R), corresponding to B22 = 1- 

6. Reduced moments i3„„ for n = 1, 2, 3 as a function of the restitution coefficient at </) = 0.2 and N = 5041. 
A similar behaviour is observed at lower densities. 

7. Reduced moments Bom and Bno as a function oi a a,t (p — 0.2 and N = 5041. 

8. Values of different coUisional averages obtained in MD simulations, as a function oi a a,t <j> — 0.05 (see 
section II B for definitions). For a < 0.5 the random rotation introduced to avoid inelastic collapse has a 



maximum deviation angle of 2.5°. The symbols 6„m are defined in Eq. (45). 
9. Distribution of (ci • C2){r) as a function of the distance between the particles at cf) = 0.05. 

10. Mean velocity-velocity correlation function at contact, as extrapolated from (ci • C2)(r) (previous figure), 
compared with the dynamic analogs (ci • C2)dyn and (ci • C2)Jj^, defined in (A. 8) and ( [A.9| ). 

11. Distribution of the impact parameter b for different a-values for cj) = 0.05. 

12. Distribution of the relative orientation of the velocities at collision (cos ■012 = Ci • C2) at a packing fraction 
of 5%. 

13. Distribution of relative velocities, Ci • C2, of colliding inelastic disks, aX (p — 0.05. 

14. To illustrate the slow reorganization of density inhomogeneities, four consecutive snapshots of the system 
are shown at a = 0.1, (p — 0.2 and N — 5000 (the full simulation box is displayed). The time interval 
between two consecutive snapshots corresponds to 50 collisions per particle. 

15. Snapshot of a typical instantaneous configuration of the system at a = 0.2, (p — 0.2 and A^ — 5000. To 
illustrate the existence of cold dense inhomogeneities, on the left (right) the particles with a less (more) 
than median kinetic energy £* are shown at real scale (i.e. the cutoff £* is chosen such that there are 
exactly half of the particles on each graph). Lengths on the x and y axes are expressed in units of the 
simulation box length. 

16. Velocity distribution of the colliding particles at = 0.05 and N = 5041: a) Original distribution; b) 
scaled velocity distributions as a function of the rescaled velocity c/ ^jT{a) for different values of a. 
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